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Abstract. The general relativistic version is developed for Robertson's 
discussion of the Poynting-Robcrtson effect that he based on special relativity and 
Newtonian gravity for point radiation sources like stars. The general relativistic 
model uses a test radiation field of photons in outward radial motion with zero 
angular momentum in the equatorial plane of the exterior Schwarzschild or Kerr 
spacetime. 



PACS number: 04.20.Cv 

1. Introduction 

For a small body like a grain of dust orbiting a star, the radiation pressure of the 
light emitted by the star in addition to the direct effect of the outward radial force 
exerts a drag force on the body's motion which causes it to fall into the star unless 
the body is so small that the radiation pressure pushes it away from the star. Called 
the Poynting-Robertson effect, it was first investigated by J.H. Poynting in 1903 pQ 
using nonrelativistic physics and Newtonian gravity and then later re-calculated in 
1937 using special relativity and Newtonian gravity by H.P. Robertson [2], who also 
stated the leading general relativistic correction to his slow motion calculation, namely 
the perihelion precession for the quasi-Newtonian orbits. These calculations were 
revisited by Wyatt and Whipple in 1950 [3] for applications to meteor orbits, making 
more explicit Robertson's calculations for slowly evolving elliptical orbits and slightly 
extending them. 

The drag force is easily naively understood as an aberration effect [4]: if the 
body is in a circular orbit, for example, the radiation pressure is radially outward 
from the star, but in the rest frame of the body, the radiation appears to be coming 
from a direction slightly towards its own direction of motion, and hence a backwards 
component of force is exerted on the body which acts as a drag force. If the drag force 
dominates the outward radial force, the body falls into the star. For the case in which 
a body is momentarily at rest, a critical luminosity similar to the Eddington limit for 
a star [5] exists at which the inward gravitational force balances the outward radiation 
force, a critical value separating radial infall from radial escape. Similarly for a body 
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initially in a circular orbit, there are two kinds of solutions: those in which the body 
spirals inward or spirals outward, depending on the strength of the radiation pressure. 

We now consider this problem in the context of a test body in orbit in a spherically 
symmetric Schwarzschild spacetime without the restriction of slow motion, and then 
in the larger context of an axially symmetric Kerr spacetime while developing the 
equations for a more general stationary axially symmetric spacetime. The finite size 
of the radiating body is ignored. Guess took this into account to extend Robertson's 
calculations in 1962 6 , and Abramowicz, Ellis and Lanza generalized this for purely 
radial motion in the exterior Schwarzschild spacetime to model jets and solar winds, 
using numerical and qualitative methods to classify the trajectories of the test body 
into seven different categories [7], introducing the key idea of a "saturation velocity" 
for which the test body remains momentarily unaccelerated due to the balancing 
of the gravitational and radiation forces. Miller and Lamb then considered the 
case of arbitrary motion around a slowly rotating star [5J [9] , after reconsidering the 
nonrotating case [TO] . 

In order to see the key features of the general relativistic case without the 
additional complications required to take into account the finite size of the radiating 
body, we return to the simple Robertson scenario and look only at how strong 
gravitational fields and arbitrary motion affect the problem. The photon flux from 
the central body is modeled by test photons in outward radial motion with respect to 
the locally nonrotating observers, namely photons with vanishing conserved angular 
momentum. While developing the basic equations for planar motion of both the 
test body and test photons in a stationary axisymmetric spacetime, we consider 
their solution only for motion in the equatorial plane of the Schwarzschild and Kerr 
spacetimes. Of course it only makes sense to consider the exterior solutions for radii 
larger than some minimum radius R outside the horizon in order to model the geometry 
outside a star (or some other physical source) of radius R producing the outflow of 
radiation. 

2. Stationary, axisymmetric and reflection-symmetric spacetimes 

Using a Boyer-Lindquist-like coordinate system {t, r, 9, <p} adapted to the spacetime 
symmetries, i.e., with d t (timelike) and d<p (spacelike, with closed coordinate lines) a 
pair of commuting Killing vectors, the metric of a stationary axisymmetric spacetime 
can be expressed by a line element of the form 

ds 2 = g tt dt 2 + 2g t<j> dtd(f) + g^dcf 2 + g rr dr 2 + g ee d9 2 , (2.1) 

where all the metric coefficients depend only on r and 9, provided that the metric 
belongs to the most interesting class of orthogonally symmetric such metrics [TT| . As 
is the case for the familiar black hole spacetimes, we further require the metric to be 
reflection-symmetric with respect to the equatorial plane 9 — 7r/2, within the allowed 
range 9 E [0, n] of this angular coordinate, where (r, <j>) then behave like flat space polar 
coordinates and will be treated as such to display orbit plots below. Our calculations 
will be limited to this plane for simplicity. 

The time coordinate lines, when timelike, are the world lines of the static 
observers. The zero angular momentum observer (ZAMO) family of fiducial observers 
has instead a 4-velocity n characterized as that normalized linear combination of the 
two given Killing vectors which is orthogonal to and future-pointing, and it is the 
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unit normal to the time coordinate hypersurfaces 

n = N- 1 (d t -N*d (j> ), (2.2) 

where N = (-g u y 1/2 and = gtcj, I g^tf, are the lapse function and only nonvanishing 
component of the shift vector held respectively. Our discussion is limited to those 
regions of spacetime where the time coordinate hypersurfaces are spacelike: g u < 0. A 
suitable orthonormal frame adapted to the ZAMOs and invariant under the symmetry 
group action is given by 

ej = n , e f = —^=d r , eg = — !=<9 e , ei = * d<f, , (2.3) 

\l9rr \I9BB \/9<t>(t> 

with dual 

J = Ndt, Lu f = y/g-^dr, J = y/g^d6, ^ = ^/g^,{d(j) + N't'dt) , (2.4) 

so that the line element (|2.ip can be also expressed in the form 

ds 2 = -N 2 dt 2 + .g 00 (d0 + N^dt) 2 + g rr dr 2 + g ee dB 2 . (2.5) 

The accelerated ZAMOs are locally nonrotating in the sense that their vorticity 
vector u>(n) vanishes, but they have a nonzero expansion tensor 9(n) whose nonzero 
components can be completely described by the shear vector 0^{n) a — 9(n) a pe^ , 
namely 

0(n) = e^<8^(n) + ^(n)®e # . (2.6) 

Since the expansion scalar 6(n) a a is zero, the expansion and shear tensors coincide. 

The nonzero ZAMO kinematical quantities (acceleration a(n) — V n n and shear 
tensor) and the conveniently defined Lie relative curvature vector [131 114] only have 
nonzero components in the r-8 2-plane of the tangent space 

a(n) = a(n) r e? + a(n) eg = d+ (lniV)ef + dsQxiN)eg , 



= W e r + 6 $ (n) B eg = -^f [d?N% + dgN% , 

k(i ic ){n) = kQ ie )(n)?e? + k (lic) (n)ge § = -[d?(ln^/g^)e f + dg(\n y /g^)e § ] . (2.7) 

Here df = ef and dg = eg. In the static limit — > 0, the shear vector 9i(n) vanishes. 

Let a pure electromagnetic radiation held be superposed as a test field on the 
gravitational background described by the metric (|2.1[) . with the energy- momentum 
tensor 

T alB = $ 2 k a k , k a k a =0, (2.8) 

where k is assumed to be tangent to an affinely parametrized outgoing null geodesic in 
the equatorial plane, i.e., fc Q V Q fc^ = with k e = 0. We will only consider photons in 
the equatorial plane which are in outward radial motion with respect to the ZAMOs, 
namely with 4-momentum 

k = E(n)[n + 0(k,n)], j>(k,n)=e?, (2.9) 

where E(n) = E/N is the relative energy of the photon and E = —k t is the conserved 
energy associated with the timelike Killing vector field and L = k^ = is the vanishing 
conserved angular momentum associated with the rotational Killing vector field, while 
z/(fc, n) defines the unit vector direction of the relative velocity. For the Schwarzschild 
case, these orbits are radial geodesies with respect to the static observers tied to the 
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coordinate system, but for the Kerr case, they are dragged azimuthally by the rotation 
of the spacetime with respect to the coordinates. 

Since k is completely determined, the coordinate dependence of the quantity $ 
then follows from the conservation equations T a ^ -p = 0, and will only depend on r 
in the equatorial plane due to the axial symmetry. From Eq. (|2.8[) using the geodesic 
condition for k, these can also be written as 

V f3 (^> 2 k /3 ) = 0, k p d p <S> 2 + fc%$ 2 = 0. (2.10) 

Integrating this equation requires an expression for the divergence of k, namely 

Ed 
k a , a =E(n)—\n[NE{n)^g^\ = - —\n[Jg^g^\ . (2.11) 

$ is then easily determined taking into account that k^dp = k r df = E/Nd?, so 

$ = {geeg^r 1 '^ , (2.12) 

as we will see explicitly below for the Schwarzschild and Kerr cases. 

Consider now a test particle of mass m moving in the equatorial plane 9 = 7r/2 
accelerated by its interaction with the radiation field, i.e., with 4- velocity 

U = ~f{U, n)[n + v{U, n)\ , v{U, n) = v r ef + v ei — vsinae? + ^cosae^ , (2.13) 

where j(U,ri) = (1 — \\v(U, n)]] 2 )^ 1 / 2 is the Lorentz factor and the abbreviated 
notation v a = v(U, n) a has been used. In a similarly abbreviated notation, v = 
\\v(U,n)\ \ and a are the magnitude of the spatial velocity v(U,n) and its polar angle 
measured clockwise from the positive cf> direction in the r-cj> tangent plane, while 
v = i)(U,n) is the associated unit vector. Note that a = corresponds to azimuthal 
motion with respect to the ZAMOs (i.e., in the <j> direction only), while a = ±w/2 
corresponds to (outward/inward) radial motion with respect to the ZAMOs. 
Using the expression (|2.2|) for n leads to the coordinate components of U 

U l = 



dt 

d7 


7 

N 1 


dr 

U r = — = 
dr 


r yv r 

\fg~vr 


de 

d7 


= 0, 


u*=** = . 

dr 


Jg~H> 



N 



(2.14) 



where r is the proper time parameter along U . Solving these for the magnitude and 
polar angle leads to 

tana 




(2.15) 

Scattering of radiation as well as the momentum-transfer cross section a (assumed 
to be a constant) of the particle is independent of the direction and frequency of the 
radiation; therefore the associated force is given by 

F(^){U) a = -aP{U) a (j Tf\ V» (2.16) 

where P(U) a p — 52 + U a Uf) projects orthogonally to U. The equation of motion of 
the particle then becomes 

ma(U) = F (nd) (U) , (2.17) 
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where m is the mass of the particle and a(U) = VjjU is its 4-acceleration. 

To examine the content of this equation, in addition to the decomposition (|2.9p 
with respect to the ZAMOs, it is convenient to decompose k with respect to U as well 

k = E(U)[U + V(k, U)} = E(n)[n + v{k, n)} . 

It then follows that 

P{U)k = E{U)V{k,U), U-k = -E{U) 

so that 

^(rad)(f/) Q = - a^{P{U) a ^) (k^l/n = a \$>E{U)\ 2 V(k, U) a . 
Therefore we hnd 

\\a(u)\\ = z[*e(u)? , <u)/\Hu)\\ = v(fc, u) , 

where a = a/m. Hereafter we will employ the simplified notation 
\\u(U,n)\\ = v, 7 ([/,n)= 7 , V{k,U) = V. 
From Eq. (|2.18p , after scalar multiplication by U, one finds the relation 
E(U) = jE(n)[l - v(U,ri) ■ v(k,n)\ 
which leads to the following expression for V 
■E(n) 



V = 



E(U) 
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n + eW)"^ ^ ~ lviyU ' ^ ' 



Introducing the polar decomposition (|2.13|) of the particle velocity Eq. 
becomes 

E(U) = jE(n)[l -vcos(a - tt/2)] = jE(n)[l - v sin a] , 
so that the frame components of V = V*n + V r e? + V*e^ are given by 

- f sin a — v ~ ~ 1 ~ i 

V =^v- ; , V = — : : r — 71/ sin a, V v = — 7^ cos a. 

l — i/ sin a 7(1 — 1/ sin a) 

Note the relations 

V* = v (V f sma + 2/ cos a), 

a{Uf = v[a{U) f sin a + a(U)*v cos a] 

which follow from the orthogonality of the pairs (V, U) and (a(U), U). 

Finally a straightforward calculation shows that the frame components of the 
4-acceleration a(U) in the equatorial plane (and hence the equations of motion) are 
given by 



2.18) 
2.19) 
2.20) 
2.21) 
2.22) 
2.23) 



2.24) 



2.25) 
2.26) 

2.27) 



7 v sin a 

„2 



dv 



a{n) r + 2i/cosa(9(n) r T + 7^— = ct[$£([/)] 2 V* , 

dr 



a{Uf 

a(U) r = 7 2 [a(n) r + fc (lic ) (n) r v 2 cos 2 a + 2v cos a 9(n) r 1 



+ 7 
a{U) 6 = , 



dv 



da 



7 since- h^cosa— 

dr dr 



a[<PE{U)] 2 V r , 



dv 



a(U)® = —71/ sinacosafcnioiM 1 * + 7 7 cos a- v since— 

dr 



da 



dr 



&[<i>E(U)} 2 V^ 



(2.28) 
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The first of these equations is a linear combination (|2.27[) of the remaining ones 
due to orthogonality of a(U) and U. Using the expression for the components of V 
one obtains 

dv sina. , „. , f . A ,„ . . 

— = \ a \ n ) + 2i/ cos a a(n) ? I H 5 — (1 — v sm a j (sin a — ) . 

dT 7 9 N^gggg^ 

dCK ""V COS QIp,.^; s\/ \r 9 7 / \ TM 

— = a(n) + 2i/cosafc>(n) ; + wkn^n) 

At v 9 

A (1 — i^sina) cos a 



dr sin a 



(2.29) 



where the positive constant A is defined by 

A = d<5>lE 2 (2.30) 
in terms of which one has 

a<£ 2 E(nf = — ■ (2.31) 

These equations admit the special solutions a = ±n/2 (radial outward/inward 
motion with respect to the ZAMOs) with 



"<"> ■ A ,|1T.| 



(2.32) 



7 y/9eeg4>4> N2 
Furthermore, if v vanishes initially the condition that it remain zero is 

; A m = a(nf . (2.33) 

This condition balancing the gravitational attraction and the radiation pressure gives 
the values r( cr it) of the radial coordinate where the particle comoves with the ZAMOs 
in a circular orbit. It will be discussed below when dealing with special cases. Recall 
that here A contains a factor l/m in cf, so that for a small enough m one can always 
achieve this force balance at a given radius. 

From the introduction it is clear that circular orbits cannot exist because of the 
drag force in the azimuthal direction. A circular orbit would have constant a = 
which requires using instead Eqs. (|2.29[) , leading to 

du A A 

-r- = m v, = ~ 1 {a(nf+2v6 $ (nyWk {iic) (ny}+ .(2.34) 

As we will see in detail below when dealing with special cases, the quantity 
Aj (y/geeg<t>4>N 2 ) is a function of r alone which in turn is a constant with respect 
to r for a circular orbit. Therefore, from the second of Eqs. (|2.34[) . we would find 
a constant value for v which is not compatible with the first of Eqs. (|2.34[) unless 
v = 0, which means U — n, i.e., the motion is comoving with the ZAMOs. The only 
surviving equation is associated with the radial component of the acceleration reducing 
the situation to the previously discussed radial force balance. In fact if one chooses 
the speed so that the radial acceleration is initially zero, then the speed will decrease, 
causing the centripetal acceleration term to decrease so that the radial acceleration 
will become negative, and orbit will decay. 
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3. Explicit spacetimes 



It is helpful to examine in sequence first the special relativistic effects involved in 
this drag force, then the addition of the gravitational field of a nonrotating source by 
considering the Schwarzschild spacetime, and then the addition of the source rotation 
finally in the context of the Kerr spacetime. The equations of motion are perhaps best 
examined in detail directly in terms of the coordinate variables t, r, <fi. The nonzero 
coordinate components of the four acceleration are 

i> dr 



a(UY 

a(uy 



d 2 t „ dt dr . . x 

d^ + 2 ^ a{n) ^ +2 -N dr 



d^r 
dr 2 

+ " 



+ (d r In y/g^) 



/dr 
Vd7 



+ 2 



\[9H> dt 
g rr dr 



dr dr 



dr dr / 



.dt d<p\ 



N 2 



N ^ + aZ ) k ^)(n)f + -f= a(n) r ( — 



dr dr / 



dt 



) Ne^ny 

2 



a(uy 



d 2 cj) yfgVr . . y 



dr 



N 



dr dr 



- 2N^ r a(n)^ _ 2^ r k {lie) (n)^ (n*£ + g) . (3.1) 

In the case of a static spacetime where iV^ = 0, 6^(n) r = 0, the above relations reduce 
to 

, TT , f d 2 t „ dt dr , , , 

° ( } = dr^ + drdr^^' 



/dr 
^dr 



fc(lie)(n)# 



iV 2 /dA 2 
-==a(«)f -T- , 



dr/ 



dr d(/> 



°(^) = dr 2 ~ 2Vft7fc ( iie)(n)f dr dr 



(3.2) 



The nonvanishing coordinate components of the radiation force (per unit mass) 
in the general case are given by 



Y 

•^(rad)^)' = 



I r -VI ^1 A 



Y 



dr 



.V" f- f^J.V 



(3.3) 



where 



, T dt , dr 

dr dr 



AX 



In the static case only the <j> component simplifies slightly due to = 0. 



3.1. Flat spacetime 

Consider first the flat spacetime case for which the metric functions are 

5« = -l, 5*0 = 0, g rr = l, gee=r 2 , g u = r 2 sin 2 



(3.4) 
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so that N = 1 and jV^ = 0. In this case, the ZAMOs are aligned with the coordinate 
time world lines and form a geodesic (a(n) = 0) and nonexpanding (9(n) = 0) 
congruence. The Lie curvature vector associated with the <f> coordinate lines only 
has a radial component with value k(\i e \(n) r = — 1/r. Moreover, radially outgoing 
(geodesic) photons in the equatorial plane have 4-momentum 

k = E{d t +d r ), (3.5) 

Thus we find E(n) = E and <f> = $ / r - 
Therefore Eqs. (I2.29[) reduce to 

dv A{1 — i/ sin a) (sin a — v) 
dr r 2 

da - Ail - sin a) cos a 

— = - i^7COsafc ( i ic )(n) H 5 . 3.6) 

dr v 

Furthermore Eqs. (|2. 14|) imply 

dr . d(j) 71/ dr 

— =7^ sin a, — = — cos a, — =rtana. (3-7) 
dr dr r dcp 

In the special case sina = ±1 of purely radial motion (+: outward, — : inward), 
one finds 

dv M±rvf dr 

- r =± , — = ±71/, (3.8) 

dr r z dr 

which leads to speeding up the outward motion and braking the inward motion as 

expected. This orbit is described by 

^ = ^£. 0.9) 
dr r^vj 

which can be easily integrated to yield a cubic relation between v(r) and r. 

Consider the behavior of the speed for purely circular initial motion sina = 
where 

dv Av _ 

J7 = --^T = - F (drag) (3-10) 

gives the magnitude of the tangential drag force per unit mass on the orbit. If one 
takes the nonrelativistic Keplerian speed v = (M/r) 1 / 2 <C 1 of a circular orbit due to 
a central mass M in Newtonian gravity then 

1/2 

%a g) = £(j) = AM <- 5/2 ■ (3-11) 

The Newtonian gravitational force per unit mass = M/r 2 only grows like the 
inverse square of the distance, so as one approaches the central mass, the drag force 
becomes more and more important compared to the gravitational free fall behavior of 
the particle initially in circular motion. Similarly the radial radiation pressure force 
per unit mass under these conditions is just A/r 2 and the ratio of the drag force to 
the radial pressure force is just namely very small. This is the approach taken 

by Robertson in describing the radiation using special relativity and gravitation using 
Newtonian theory. 

Since we use geometrical units c = 1 = G, to compare these formulas with the 
literature, one must restore the standard units. If m is the mass of the particle under 
consideration, then the tangential drag force is actually 

W 

m -F(drag) = -gv (3.12) 
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where v = (GM /r) 1 / 2 is the Keplerian speed and W = cL a / {Airr 2 ) is the power of 
the incoming radiation at the particle position expressed in terms of the luminosity 
C of the source and the cross-section a of the particle. Introducing the Eddington 
luminosity 

4nGMcm v , R / M \ , 

£ E dd = P - = 1.3 x 10 38 — erg/s , (3.13) 

where m p is the proton mass, ctt is the Thompson cross-section and Mq is the mass 
of the sun, the drag force can be rewritten 

£_\ /mp\ fo_\ {GM) 3 / 2 
*Edd / V m / I 0"T , 



mf «= m z^K^J UrTv^- (3 - 14) 



This identifies our radiation constant as 
A = GM ' £ ' / 7 



(^) • (3.15) 
v m J 



. ^Edd / V °t / \ m 

If we consider a hydrogen atom as representative of the particle, then its mass is m p 
and <tt = c leading to the ratio 

^ - (a;) (3 - 16) 

which is small for sources of astrophysical interest. 
3.2. Schwarzschild spacetime 

The Schwarzschild spacetime is characterized by the metric functions 

3tt = -iV 2 , g t<f> = 0, g rr = l/N 2 , g eg = r 2 , g^ = r 2 sin 2 9 , (3-17) 

where the lapse function is 



2M M , 

N=\\ 1- — , (3.18) 

V r r 

in which the approximate expression represents the asymptotic value at r — > oo and to 
first order in M. In this case as well the ZAMOs are aligned with the coordinate time 
world lines; however, they form an accelerated (a(n) r = M/ [r\Jr 2 — 2Mr) ~ M/r 2 ) 
and expansionfree (6{n) = 0) congruence. The Lie curvature of the <$> loops has only 
a radial component with value fc(iie)(»i) r = —N/r ~ — 1/r + M/r 2 . Radially outgoing 
(geodesic) photons on the equatorial plane have 4-momentum 

-l 



E 



d t + d r 



E(n)[n + e f }. (3.19) 



Eqs. ([2~29f reduce to 



dv iVsina 2 A 

— = v K + -~^{l-V8ina){ama-v), 

dr r7 r^SS 1 

(v 2 — v K ) + 9 (1 — 2/ sin a) cos a , 



da Njcosa A 



dr 

— =7^iVsina, (3.20) 
dr 

where we have used the relation a(n) r = — fc(ii )("-) r ^X' 
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is the Keplerian speed associated with circular geodesies. If one is only interested in 
the spatial orbit of the particle, one can choose r or as the parameter along its path, 
re-expressing the above equations using the chain rule. 

For the case v = of a particle at rest with respect to the ZAMOs, these equations 
reduce to the single condition Eq. (|2.33|) representing the balancing of the gravitational 
attraction and the radiation pressure at constant r and 0, namely 

A / 2M\ 1/2 2M 

M = { 1 -—) "> r = r ^ = 1 - A 2 /M 2 ■ (3 ' 22) 

This is exactly the zero speed singular point discussed for radial motion taking into 
account the finite radius of the photon source by Abramowicz et al [7] (which therefore 
introduces a complicating geometric factor into the problem incorporating the effective 
solid angle of photons which arrive from the entire surface of the photon source) and 
makes the general relativistic case quite different from the Robertson limit where the 
balance condition is independent of radius. Instead for the general relativistic case in 
the present context, if A/M > 1, the net effect is always expulsion from the central 
object as in the Robertson limit, but if A/M < 1, there is always exactly one radius 
rb S (2M, oo) at which force balance is achieved for zero speed, and this point in 
the system of differential equations behaves like an attractor for some subspace of the 
initial data. If A/M <C 1, then this radius is near r = 2M and hence not important for 
a real star. In fact in the Newtonian limit this condition reduces to simply A/M = 1 
which is independent of the radius. 

This equilibrium condition can also be re-expressed in terms of the Eddington 
luminosity as above using its general relativistic generalization 7J 

£ E ddGR = £ E dd/(l - 2M/R) 1 / 2 , (3.23) 
where R is the radius of the emitting body. Then the equilibrium condition is 



.4 



( C \ f l-2M/r \ 1/2 fa_\ /m^ 
V^EddGiJ \1-2M/RJ \ar)\m) 



GM V^EddGR./ Vl '' ~ M '' 

so that when evaluated for a proton at the surface of the star where the Thompson 
cross-section is relevant, one has a balance of the outward radiation pressure and the 
inward gravitational force. 

The system of differential equations (|3.20p requires numerical solution, easily 
done with a computer algebra system (Maple was used here). One can study the 
initial behavior for a solution which starts out close to a circular geodesic at a fixed 
radius rg and constant speed vk to see how the orbit begins to decay for the more 
interesting case A/M < 1. This is done in appendix B. 

Figs. 1 and 2 show some typical solution curves starting from initially circular 
initial data either inside (Fig. 1) or outside (Fig. 2) the critical radius at which a 
particle at rest with respect to the ZAMOs (which in turn are at rest with respect 
to the coordinate system) remains at rest. For comparison, geodesies with the same 
initial data are shown in gray. Unless the particle has a sufficiently high initial speed 
that it can escape to infinity, it is forced to migrate to this critical radius where it 
comes to rest. Fig. 2 shows two final cases where the initial speed is sufficient for either 
the solution curve or the geodesic to escape. If the initial radius is taken far outside 
the critical radius, the circular orbit gradually decays as in the Robertson limit until 
it finishes at the critical radius. 
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Figure 1. The orbit of the particle in the Schwarzschild spacetime with M = 1, 
A/M = 0.6, r (crit ) = 3.125A/, v K = 0.7071. The inner circle is the horizon 
r = 2M, while the outer circle is at the critical radius which is inside the 
initial data position. Initial conditions have (r(0), </>(0), «(0)) = (AM, 0, 0) and 
^(0) = 0.2,0.5,0.8. The corresponding geodesies A/M = are in gray. 




Figure 2. The orbit of the particle in the Schwarzschild spacetime with 
M = 1, A/M = 0.8, r (crit) = 5.5A/, u K = 0.7071. The inner circle is the 
horizon r = 2M, while the outer circle is at the critical radius which is outside 
the initial data position. Initial conditions have (r(0), </>(0), »(0)) = (AM, 0, 0) 
and for the left figure f(0) = 0.2, 0.3, . . . , 0.7 < vk while for the right figure 
0.71, 0.72, . . . , 0.75 > vk- The corresponding geodesies A/M = are in gray. 
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Figure 3. The orbit of the particle in the Schwarzschild spacetime with M = 1, 
A/M = 0.01, r (crit) « 2Af and initial conditions (r(0), </>(0), af(0)) = (10iW, 0,0) 
with the circular geodesic speed u(0) = vjf = 0.3536. The circular geodesic (gray) 
is shown with the in-spiraling orbit (black) with the same initial conditions. 



3.3. Kerr spacetime 

In the equatorial plane of the Kerr metric, the metric is 

/ 2M\ 2aM r 3 + a 2 r + 2a 2 M 

9tt = — 1 , 9t<t> = , 94>4> — , 

\ r J r r 

r 2 

9rr = , gee =r 2 , (3.25) 
so that 



N= / rA \-— N* - 2aM 2aM (3 2fi^ 

V r 3 + a 2 r + 2a 2 M ~ r ' r 3 + a 2 r + 2a 2 M ~ r 3 ' 1 ' J 

where A = r 2 + a 2 — 2Af r and the approximate expressions represent their asymptotic 
values (r — > oo) to first order in M. The ZAMOs are timelike outside the horizon 
r+ = M + V M 2 — a 2 . The nonvanishing components of the ZAMO kinematical fields 
are 

M[(r 2 + a 2 ) 2 - Aa 2 Mr] M 



a(n 



r 2 VA(r 3 + a 2 r + 2a 2 M) 

{nf = 



aM(3r 2 + a 2 ) 3aM 



r 2 (r 3 + a 2 r + 2a 2 M) r 3 ' 

. . f \/A> 3 - a 2 M) 1 M 

v i T 2.(ji _|_ a ^ r _|_ 2a z M ) r r z 

Circular geodesies correspond to orbits 

a 2 =F 2a\f~Mr + r 2 

lT ± = 7± (n + , ± e,), I/± = ___ I (3.28) 

and the following relation between i/± and the ZAMO kinematical fields hold 

a(n) f = fe (lie) (n)V+^_ , -2^(n)' = fc (lio) (nf + (3.29) 
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Figure 4. The orbit of the particle in the Kerr spacetime with M = 1, a = 0.5 
(left figure), a = -0.5 (right figure), A/M = 0.6, r (crit) = 3.154M. The inner 
circle is the horizon r = 1.866M, while the outer circle is at the critical radius 
which is inside the initial data position, initial conditions have (r(0), (f)(0), ct(0)) = 
(4M, 0, 0) and u(0) = 0.2,0.5,0.8. The corresponding geodesies A/M = are in 
gray. The bound orbits end up co-rotating with the hole at the horizon (geodesies) 
or at the critical radius (accelerated). 



Equations (12.29)) can then be rewritten as 



du 
dr 

do 
d7 
dr 
d7 



sin a 



7 

7 cos afyii, 



v+v- — v cos a (v+ + !/_)] + 



v cos a (y+ 



A{1 — v sin a) (sin a — v) 

2 All — vsma) cos a 
hi/ + ; AT2 



sm a 



(3.30) 



For the case v = of a particle remaining at rest with respect to the ZAMOs, 
these reduce to the single radial force balance condition 



A 
M 



[(r 2 



,2\2 



4a 2 Mr]VA 



! (S^) 3/2 



(3.31) 



which cannot be solved explicitly for r. However, the right hand side of this equation 
takes values between at the horizon r = r + and 1 when as r — > oo so a critical radius 
rv cr it) always exists for which this is satisfied for any proper fractional value of A/M. 
If A/M > 1 of course no static solutions exist. This balance condition can also be 
used to generalize the expression for the Eddington luminosity to this case and rewrite 
the balance condition in terms of it as done for the Schwarzschild case. 

Fig. 3 shows some typical solution curves starting from initially circular initial 
data outside the critical radius at which a particle is at rest with respect to the ZAMOs 
in the case of positive rotation parameter a — 0.5 corresponding to counterclockwise 
dragging (left) and in the case of negative rotation parameter a = —0.5 corresponding 
to clockwise dragging (right). The same qualitative orbit behavior results as in the 
Schwarzschild spacetime but here the ZAMOs themselves are dragged along by the 
spacetime geometry, so that the actual curves in the plane are severely distorted by 




Figure 5. The orbit of the particle in the Kerr spacetime with M = 1, 
a = 0.5 (left figure), a = -0.5 (right figure), A/M = 0.8, r (crit) = 5.551M. 
The inner circle is the horizon r = 1.866A/, while the outer circle is at the 
critical radius which is outside the initial data position. Initial conditions have 
(r(0),</>(0),a(0)) = (4Af,0, 0) and u(0) = 0.2,0.5,0.8 for the left figure, while in 
the right figure v{0) = 0.2,0.5,0.8,0.847 for both the accelerated and geodesic 
curves and then finally f(0) = 0.9 for the accelerated curve and i/(0) = 0.848 
for the geodesic, both of which escape to infinity. The corresponding geodesies 
A/M = are in gray. The bound orbits end up co-rotating with the hole at the 
horizon (geodesies) or at the critical radius (accelerated). 

this effect. Unless the initial velocity is sufficiently high, the solution curves end up at 
corotating with the geometry at the critical radius. Fig. 4 shows the same situation 
for the case in which initial data is taken inside the critical radius. Initial data is 
taken at rather small radii not far from the horizon only to show the effects of the 
radiation pressure in their most exaggerated form, and is not intended to model actual 
astrophysical scenarios. 

4. Concluding remarks 

It is remarkable that in the seventy years since Robertson studied this effect, there 
seems to be no published article which reconsiders it within the full theory of general 
relativity. It is a clean geometric model that is a simple extension of geodesic motion 
and with today's computer algebra systems, one can plot orbits accurately within a 
matter of seconds to explore its consequences numerically. Because Robertson was 
himself repeating earlier nonrelativistic calculations, his own article omitted many 
details which we have provided here in the framework of stationary axisymmetric 
spacetimes, made explicit for the Schwarzschild and Kerr spacetimes. The model 
we developed here, while still a toy model, can provide interesting applications to 
astrophysical problems, especially in relation to accreting black holes and neutron 
stars. The inclusion of effects such as non-equatorial particle orbits and photon fields 
endowed with angular momentum will be especially relevant in this context. 
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Appendix A. Weak field, slow motion, small drag limit 

Robertson developed the relativistic equations appropriate to describe this problem 
in special relativity and then looked at the Newtonian approximation, taking into 
account by hand the perihelion precession due to general relativity. His results can be 
obtained from the Schwarzschild case in the weak field, slow motion approximation 
for a small drag coefficient. Letting / = df/dt, and introducing the approximation 



t 

one obtains 



M 1 A 1 
< 1, — < 1 , — < 1, 
r r 



The azimuthal equation leads to the constant of the motion 

h = r 2 (j) + A<f> 

so that H = r 2 (f> = h — A<j). Backsubstituting this into the radial equation, introducing 
the dimcnsionless reciprocal variable u = M/r, and changing the independent variable 
from t to <p using d(f>/dt — H/r 2 , one finds 

d 2 u M(M -A) Adu , 

W +U= H 2 -Hd^ + 3U ■ 

The first term on the right hand side represents the Newtonian gravitational force 
slightly reduced by the radiation pressure which leads to elliptical orbits, the last 
term the general relativistic term responsible for the perihelion precession of those 
elliptical orbits, and the middle term the drag effect responsible for the decay of the 
orbit. 

In the absence of the precession term this reduces to 
d 2 u A du _ M(M - A) 

d^ + (h - A<fi) M> + U ~ (h- A(f>) 2 ' ( 

which is further simplified as noted by Robertson by introducing the new variable x 
such that 4> = —x + h/A 

d 2 u 1 du q M(M -A) fl n . 

which admits the following solution in the notation of Maple 

u = x[Ci BessclJ (1, x) + C 2 BesselY (1, a;)] + q xLommclSi (-2, 1, x) , (A. 3) 

with the consequence 

du 

— = x[CiBesselJ(0,a;) + C 2 BesselY (0,x)} - 2qa;LommelSi (-3,0, x) . (A.4) 

This exact solution u(x) of the approximate equations can be used to express a(x) 
and in turn v(x). In fact, using the approximate relations 

dr d(j) v 

— = vsma, — — — — cosa, (A. 5) 

At At r 1 
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one gets 

dr du 

— = rtanct — > cot a = it— , (A. 6) 

d0 dx 

where du/dx is given by Eq. (|A.4[) . 

However, one can see how the approximately elliptical orbits evolve including the 
relativistic precession term by using the further approximation of small drag A/h <C 1, 
so that 

1 If A 



Then radial equation linearized in this ratio leads to 

d 2 u M(M — A) f 24 A 2 A du , . « 

Solving this in the absence of the A/h and precession correction terms leads to a 
Newtonian elliptical orbit 
M 

U( Q ) = — (l + ecos0), (A. 9) 

where the orbital parameters correspond to an adjusted central force M — > M — A 

P = a(l-e*) = v Jt Ty (A.10) 

Excluding the single nonlinear precession term and using an osculating ellipse 
approach to investigate the secular perturbations of these elliptical orbits to first order 
in A/h and to first order in u <C 1, one looks for approximate next order solutions 

U{4>) = U(Q){4>) +M(1)(0) 

and finds 

M ( 1 N 
= I 20 — -e0cos</> 

Recombining the various terms to first order in A/h leads to 

M 

u(<t>) = -rr?{l + e(cj>) cos 0-5(0)). 

where 

p(0) =p(l- 2^ , e(0) =e(l- , (5(0) = , 

which correspond to the following differential equations satisfied by the orbital 
elements for their logarithmic rates of change and in turn that of a = p/(l — e 2 ) 
using logarithmic differentiation 

dlnp A dine 5 A dlna 2 + 3e 2 A 

d4~ = ~ 7i' d(f) ^"Yh 1 d<\> = ~ 1 - e 2 ~h ■ 

To convert the derivative into a t derivative in these equations, one can take the 
average related rate over one period of the zeroth order solution 

dt _ r 2 _ M 2 dt (0) _ M 2 p 2 [o(l - e 2 )] 2 



H Hu 2 d(j) hu 2 0) h(l + ecos0) 2 h{l + > 
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by calculating 

dt\_ At\ A ^ _ [q(l - e 2 )] 2 /** dcj> _ a 2 (l - e 2 ) 1 / 2 



2tt 2tt/i J (1 

from which Robertson's results follow in the form similar to that given by Wyatt and 
Whipple [3] for the secular changes in these orbital parameters 

dlnp A dine _ 5^4 dlna 2 + 3e 2 A 

eft a? ' eft = ~2^ ' eft = ~ 1 - e 2 a 2 " ' 

The quadratic precession term can then be handled again by considering a second 
linearized term U( 2 ) (</>) in the solution, following the discussion of Ohanian and Ruffini 

u{4>) = "(O)(0) + U(l)(4>) + «(2)(0) ■ 

Expanding the 3m 2 term about the zeroth order solution 3u 2 -j in the original second 
order differential equation (|A.8[) . one first neglects the small uru and U( 2 ) correction 
terms in the square compared to 3m 2 ^ , leading to 

d 2 U(2) 2 

~W + U(2) = 3 » 

M 2 

= 3 — ^-(1 + 2e cos + e 2 cos 2 0) 

M 2 / 1 1 \ 

= 3— 1 + -e 2 + 2e cos + -e 2 cos 20 . (A.ll) 

p V 2 2 / 

The cos 20 term on the right hand side leads to small oscillations in the solution for 
U( 2 ), and the constant term a slight constant shift, but the cos0 term is in resonance 
with the left hand side oscillator and will dominate both over the long term, leading 
to the secular part of the solution 

M 2 e 

U(2)(<£) =3—5-081110. (A.12) 

Reconstructing u to first order in this small quantity then leads to the same form as 
above but with the additional precession term 

5U) =3—0, (A.13) 
P 

which agrees with Robertson's result when A -C M. 



Appendix B. Perturbing initially circular orbits in Schwarzschild 
spacetime 

To consider an orbit in the Schwarzschild spacetime close to a circular geodesic 

a (r)=0, v{t) = v k , r(r) = r . (B.l) 

at a fixed radius ro and constant speed vk for the case of a small radiation pressure 
term whose effect will be small at least initially, we linearize Eqs. (I3.20[) with respect 
to the dimensionless parameter e = A/M which is assumed to be much less than 1 by 
expanding those equations as follows 

a (r) = eai(r) , v(t) = v K + ez/i(r) , r(r) = r + e n(r) . (B.2) 
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Inserting these quantities into Eqs. (|3.20[) and linearizing them in e we find 

— = ax , 

dr 7k r 

dai 2( K j K VkCkIK , v K 
dr vk M ro 

"j~ = IkCkTooii , (B.3) 
dr 

where (k = (M/r^) 1 / 2 is the angular velocity associated with the circular geodesic. 
We solve these with initial conditions representing the initial circular geodesic motion 

n(0) = 0, ai(0) = 0, i*(0) = 0, (B.4) 
leading immediately to 

£(0) = 0, ^(0) = ^, = (B.5) 

dr dr tq dr tq 

Differentiating the second of Eqs. (|B.3|) and using the other two gives a decoupled 
equation for a± 

d 2 ai 
dr 2 

where 



n A ai =b, (B.6) 



Note that f2 is a real frequency when ro > 6M and coincides with the so called 
"epicyclic frequency," rescaled to correspond to a proper time parametrization of the 
orbit [16] , The solution for a\ is given by 

«i( r ) = TS~ I n sin ( 0r ) + 27Jr^C/c(co8(nr) - 1)] , (B.8) 



Introducing the notation 



S(j)= f T ai(r')dT' 
Jo 



n 2 r Q 

Then one immediately gets 



(1 - cos(ftr)) + 2v K "i K C K ( — sin(fir) - r 



(B.9) 



ri(r) = 7i^oS(r), ^i(r) = — — — — — iS'(r) — ^-r . (B.10) 

Introducing the geodesic proper time orbital frequency fi( or b) = Jk^k/to, the 
oscillation frequency for the perturbation can be expressed as 



"(orb) V T 

The final nonoscillatory term in S(r) proportional to r leads to the decrease of the 
radius causing the orbit to decay. 
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